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Abstract 

We analyze the microscopic dynamics and transport properties of a gas of thin hard rods. Based 
on the collision rules for hard needles we derive a hydrodynamic equation that determines the 
coupled translational and rotational dynamics of a tagged thin rod in an ensemble of identical 
rods. Specifically, based on a Pseudo-Liouville operator for binary collisions between rods, the 
Mori-Zwanzig projection formalism is used to derive a continued fraction representation for the 
correlation function of the tagged particle's density, specifying its position and orientation. Trun- 
cation of the continued fraction gives rise to a generalised Enskog equation, which can be compared 
to the phenomenological Perrin equation for anisotropic diffusion. Only for sufficiently large den- 
sity do we observe anisotropic diffusion, as indicated by an anisotropic mean square displacement, 
growing linearly with time. For lower densities, the Perrin equation is shown to be an insufficient 
hydrodynamic description for hard needles interacting via binary collisions. We compare our re- 
sults to simulations and find excellent quantitative agreement for low densities and qualtitative 
agreement for higher densities. 
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I. INTRODUCTION 



Theoretical work on the transport behavior of non-spherical molecules was pioneered by 
the work of Debye^ and Perrin 2 ^. Brownian motion of a rotational ellipsoid in a viscous 
solvent is described in terms of a partial differential equation coupling rotational and trans- 
lational diffusion (see also Ref. 0). The central quantity is p(r, u, t; r', u', 0), the probability 
to find a molecule with center of mass (CMS) position r and orientation u at time t, given 
that its CMS position is r' and its orientation is u' at time t = 0. Its dynamical evolution 
is governed by the following translational-rotational diffusion equation, first introduced by 

p err j p 2 J 



d t p(r, u, t; r', u', 0) = [-D R L 2 + D\\ (V r • u) 2 + D ± (V 2 - (V r • u) 2 )] p(r, u, t; r', u', 0). (1) 

Here u denotes the orientation of the symmetry axis of the molecule. The operator L 2 is 
the familiar angular momentum operator. It describes rotational diffusion, Dr being the 
rotational diffusion constant. The usual Laplacian has been decomposed into directions par- 
allel and perpendicular to the particle orientation, with parallel and perpendicular diffusion 
constants D\\ and D± respectively. 

Transport coefficients may be calculated using the Kirkwood-Riseman theory^ for rod-like 
molecules in solution. Broersma£ computes the rotational and isotropic translational diffu- 
sion coefficients approximately for long cylinders with no stick boundary conditions. Rods 
with capped hemispheres at the ends have also been considered^. Alternatively the rods 
are composed of spherical hydrodynamic subunits ^ 10 ' 11 with each sphere interacting with 
the other spheres via hydrodynamic interactions. The resulting transport coefficients obey 
the Stokes- Einstein relation, i.e. both the rotational and the translational diffusion constant 
are inversely proportional to the solvent viscosity. Transport properties of nonspherical 
molecules have also been computed using kinetic theor y 12 ' 1 ?' 14 ' 1 ^ 

In polymer physics, the transport behavior of rod-like macromolecules has been investi- 
gated where anisotropic translational diffusion and the coupling of rotational and transla- 
tional diffusion is most pronounced. The coupling of rotational and translational diffusion 
can be detected experimentally by depolarized light scattering. Dynamic light scattering 
from cylindrically symmetric rigid particles has been discussed by Aragon and Pecor a 1 ^ 1 ^ 1 ^ . 
If only autocorrelations are included, the intensity of scattered light is proportional to the 
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intermediate scattering function 

1 N 

i=l 

3 

with a w = n T a % n f- 

Here k is the wave vector and a^l is the polarizability tensor of particle i and n mc (n fin ) 
denotes the polarization vector of the incident (scattered) light. Aragon and Pecora have 
solved the Perrin equation, allowing for a comparison of the theoretical predictions of the 
Perrin equation and the measured intensities. If the polarization of the incident and scattered 
light are orthogonal, then the scattering intensity is proportional to 

^ N ■ 

^(k,t)~^E((^(°)-^(*)) 2eik(rIW " ri(0)) ) ( 2 ) 

i=l 

If the average over orientation and position is factorized, the intermediate scattering function 
decay with a rate T = 6Dji + D iso k 2 . The factorizing assumption, however, is valid only if the 
respective decay rates for orientational und positional degrees of freedom differ sufficiently 
This is usually the case as long as > D iso k 2 , which is guaranteed for small wave vectors. 
Hence translational and rotational diffusion coefficients can and have been obtained from 
polarized light scattering experiment o 19 i 20 i 21 , however the experiments are difficult due to 
low intensities of the scattered light^ 2 -. Since rotational and lateral motions of a rod-like 
molecule are usually studied in light scattering experiments for kL > 1, where L is the rod 
length^ 2 ., the factorizing assumption is bound to break down and coupling of translational 
and rotational diffusion will become important. 

Even though the Perrin equation was originally formulated for the Brownian motion of 
a rodlike macromolecule in solution, it is usually thought to be applicable to the motion of 
a tagged rod in a fluid of otherwise identical rods. These systems have been simulated in 
the limit that the diameter of the rods goes to zero. Frenkel and Maguire 2 ^ suggested an 
algorithm for the dynamics of hard needles and furthermore computed the autocorrelation 
of the linear velocity and the angular velocity in the Enskog approximation. The latter 
breaks down when the density po is sufficiently high, such that the tagged rod cannot rotate 
unless other rods move out of the way. This so called semidilute transition density was 
estimated to be^ p = yL 3 ~ 70. Here N denotes the total number of rods of length L 
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confined to a volume V. For larger densities the rotational diffusion constant Dr ~ p$ 2 , 
as first suggested by Doi and Edwards with help of a scaling argument^. They argue that 
rotational diffusion of a rod can be described by a random walk, so that the rotational 
diffusion constant is given by the product of the jump frequency, u, and the average step 
size squared: Dr = viajV) 2 ■ In the semidilute regime the jump frequency is v ~ D\\/L 2 
and the step size is (a/L) ~ 1/poj leading to the above scaling of Dr with density. A similar 
scaling argument has been formulated for the transverse diffusion coefficient. The inverse 
of the jump frequency v can also be interpreted as the disentanglement time, or within the 
reptation model as the time needed for the breakup of the tube, confining the test rod. 
On time scales of order \jv the rod can diffuse in the transverse direction by an amount a, 
giving rise to D± ~ a 2 is ~ p$ 2 . Such a scaling has been derived by Szamel and Schweize r 26 i 27 
within a dynamic mean field theory for self and tracer diffusion. They neglect rotational 
motion of the rods and close the hierarchy of dynamical distribution functions on the two 
particle level. 

In this work, we discuss the microscopic dynamics of a gas of infinitely thin needles as 
a minimal model for rod-like macromolecules. An approximate kinetic theory is derived 
which allows us to determine the time delayed autocorrelation function of a tagged par- 
ticle's density with a given orientation u. We derive a hydrodynamic equation for this 
autocorrelation function from the microscopic collisional dynamics and examine the validity 
of the phenomenological Perrin equation in the limit of small wave vectors kL <^ 1. In fact, 
differences with the Perrin theory are found. In particular, meaningful anisotropic diffusion 
constants D\\ and D± can only be obtained in the limit of high densities. In contrast to the 
Perrin theory we include the free motion of the particles in the microscopic theory, giving 
rise to ballistic behavior at short times. In addition to smooth needles, we consider rough 
needles^, which allow for a change of the relative tangential velocity in a collision. The 
results, however, are qualitatively the same as for smooth needles. 

The outline of the paper is as follows. In Sec. |H] we analyse the solutions of the Perrin 
equation following Aragon and Pecora. In particular we discuss anisotropic diffusion with 
help of the mean square displacement parallel and pependicular to the rod's orientation. In 
Sec. IIHI we analyse binary collisions of two hard needles and introduce a Pseudo Liouville 
operator for the time evolution of a gas of hard needles. Subsequently (Sec. HVjl we apply the 
projection operator formalism of Mori and Zwanzig to the correlation function of interest, 
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which is thereby represented as a continued fraction. The formalism is applied to the density 
of a tagged particle with given orientation, i.e. p(r,u), in Sec. El In Sec. IV CI we describe 
numerical simulations. The results of a truncation of the continued fraction are presented 
and compared with the simulation results in Sec. IV1I We conclude with a discussion and 
outlook in Sec. IV11I 



II. ANALYSIS OF THE PERRIN EQUATION 

In this section we review the predictions of the Perrin equation, Eq. (HI), as discussed in 

nnn 

Refs. I16II17U18I . with particular emphasis on anisotropic diffusion. We do not present any 
new results, instead we want to familiarize the reader with this classic phenomenological 
theory. Later on, we will compare this theory with our kinetic theory, which is based on a 
microscopic dynamic model. 

The first step in the solution of the Perrin equation is a Fourier transform with respect 
to the spatial coordinates 

p(k, u, u'; t) = J d 3 r e ikr p(r, u, t; r', u', 0) (3) 

where we have assumed spatial homogeneity. The Fourier transformed function is the solu- 
tion of 

d t p(k, u, u', t) = i-D R L 2 - k 2 Dj_ - pi, - Dj_)(k • u) 2 ]p(k, u, u', t). (4) 

It is convenient to choose a coordinate frame k = kz with k — |k|, so that the unit vector u 
is specified by a polar angle <p and r\ — k- u/(k). In terms of 4> and rj the angular momentum 
operator is given explicitly by: 

~ 2 d 2 d Id 2 

drj^ ^ ^ dr/ 1 — rf- dcfi 2 ^ ^ 

Provided 7 2 := (D\\ — Dj_)k 2 / Dr > 0, the right hand side of Eq. (J3J) can be identified with 
the anisotropic Laplacian, whose eigenf unctions, S mi i('y,ri)e' im ^, are known and expressed in 
terms of the prolate spheroidal wavefunctions Si m 29 '* 30 . The latter solve the equation 
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-^(1 - r, 2 )^ + ^ + 7 V ) S lm (%V) = \ lm (l 2 )S lm (l,v) (6) 
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with eigenvalues 

^ - + + Iw 1 ' + ohl) (7) 

We are interested in diffusion and hence consider only long wavelength phenomena, so that 
it suffices to keep terms up to and including 0(j 2 ). 

The full solution of the Perrin equation can then be written in terms of the eigenf unctions 

1 oo 

p(k, u, u', t) = — 2 £ E S ^ il)SUl, n'y^e-^'e-^ 2 * (8) 

m=0 l=m 

such that the eigenvalues determine the relaxation rates according to T/ m (7 2 ) = D±k 2 + 

The anisotropic diffusive behaviour can be analysed and has been discussed in the lit- 
erature in terms of the mean square displacement. We define the displacement Ar(t) := 
ro(t) — ro(0) and compute its projection on the initial orientation Ar\\(t) := Ar(i) ■ Uo(0) 
as well as perpendicular to it Ar±(t) := Ar(t) — uo(0)Ar(£) • Uo(0). The Perrin equation 
predicts 

((Ar„(t)) 2 ) = 2D iso t + 2(D l~° ±] (1 " e- 6 ^) (9) 

((Ar ± (t)) 2 ) = 4Asot - 2{D 1 Dr D±) (1 - e~ 6 ^) • (10) 

Analysing the short time behaviour one finds anisotropic diffusion, ((Ar||(t)) 2 ) ~ 2D\\t 
and ((Ar±(t)) 2 ) ~ 4D±t, such that in general the diffusive motion parallel to the initial 
orientation is faster than perpendicular to it. The asymptotic long time behaviour is again 
diffusive with, however, an isotropic diffusion constant D iso = (Du +2D±)/3. The crossover 
is determined by the rotational diffusion constant Dr, such that for QD^t S> 1 the initial 
orientation has been forgotten and diffusion is completely isotropic. 

The above discussion already indicates that it might be difficult to detect anisotropic 
diffusion, since it is not the true asymptotic long time behaviour and, on the other hand, 
we expect ballistic motion for short times which is not accounted for in the Perrin equa- 
tion. Hence the anisotropic diffusion on intermediate timescales will be superimposed by a 
crossover to ballistic motion on short time scales and isotropic diffusion on long time scales. 
Only if the density is sufficiently high, so that the time for rotation of the molecule is large, 
do we expect to see a clear separation of the two diffusive regimes. 
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Several attempts have been made to detect anisotropic diffusion from velocity autocor- 
relations. For isotropic diffusion there is a simple relation, connecting the asymptotic be- 
haviour of the mean square displacement with the velocity autocorrelation, integrated over 



i)) 2 ) is not simply related to the autocorrelation of 



time. For anisotropic diffusion ((Ar 
v Jt) = v(t) • u(0), as noted in Ref. \2A 



III. DYNAMICS OF HARD NEEDLES 
A. Collision rules 

The phase space of needles is specified by the needles' translational and angular velocities 
Vj, u)i and their center of mass positions rj and orientations Uj, with % — . . . N. The needles 
have length L, mass m and moment of inertia /. Given two needles, labeled and 1, their 
orientations define a plane Eq\ whose normal vector is given by (see Fig. Q): 

u ± = 1 r (11) 

|u X Ui| 

The three vectors u , and 

U X - (U • Ui)llo 



Ur 



1 - (u • Ui 



|2 



(12) 



define a local orthonormal frame which is convenient for our calculations. For there is an 
analoguous expression if one permutes indices and 1 in the last equation. Let r i = r — ri 
be the CMS distance between two needles, which may be decomposed into a component 
perpendicular r$i = (r i • uj_)uj_ and parallel = s iu — Si U! to -E^i— . The parameters 



are determind from 



Note that uj~ = — yl — (uq ■ Ui) 2 uq + (uq ■ Ui)uq 



% = ~T^—T\ (13) 



32, with the coeffi- 



Let us recall the general collision rules for hard needles (see e.g. Ref. 
cient of the normal restitution e = 1 and the coefficient of tangential restitution (3 = ±1). 
The relative velocity of the contact points is given by V = vo — vi + soi u o ~ sioiii 
and can be decomposed with respect to the orthonormal frame (ui,ujj-,u ) according to 
V = xuj_ + y\iQ + zu . The normal and tangential component of the relative velocity after 



7 



the collision (primed quantities) are given in terms of the respective components before the 
collision as follows: 



(ul • V) = -(ul • V) (14) 
(u_l x V) = -P(u ± x V) (15) 

Here (3 = —1 for smooth needles, and — 1 for totally rough needles, being the only values 
allowed by energy conservation. In terms of the particle's translational and angular velocities 
the collision rules are the following: 

Po = Po + Ap 
Pi = Pi - Ap 

%1 A ( 16 ) 

UJ = LJ + — u x Ap 

/ S 10 A 

(jJ x = U>i — 7~ Ul x P" 

B. Perfectly smooth needles: (3 = — 1 

In this case, the translational momentum transfer occurs only perpendicular to the plane 
Eq\, and is given by Ap = au±_ with 

a = —m -. — — — v (17) 

l + g(i + ^o)) 

The temporal change of orientation u is changed by a collision due to translational mo- 
mentum transfer, / Au = SoiAp. 

C. Perfectly rough needles: (3 = 1 

Rough needles also allow for a momentum transfer in the directions uo and u^: 

Ap = au ± + 7iu + 72^. (18) 
The coefficients 71 and 72 are given by: 

7l = l±l B V- C * (19) 

71 2 AC-B^ 1 } 

l+(3Bz-Ay 

72 = —AC^& (20) 
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with 



A = - + |^(l-( Uo -u 1 ) 2 ) (21) 
m 21 

B = -^(uo.mj^l-fuo.m) 2 (22) 

and the variables y, z as defined above. The change of angular velocity is given by the 
expression: 

Au = ^(au ± + 72 u^). (24) 



D. Liouville operator for hard needles 

The time evolution of any function A(T,t) on phase space T = {r^, iij, Vj, uJi} is deter- 
mined by the pseudo-Liouville operator £ + 13 i 14 

A(T, t) = exp(iC + t)A(T, 0), t > (25) 

The pseudo-Liouville operator £ + consists of two parts, C + = C° + C + . The first part results 
from the Poisson bracket with the Hamiltonian of free translational motion and rotation, 
i£° = {Ti, . . .}. For N needles the Hamiltonian expressed in terms of generalized canonical 
momenta is given as follows: 

The second part of the pseudo-Liouville-operator, describing binary collisions, reads 



H = £ ( h< + hi + ^ji^pI ) ( 26 ) 



= \ E sKl 6 f-^l^i) 6 (i - * (f - M ) 5 d r i' - 0+ ) K - 1) 



2 . 



(27) 

Here #(x) is the Heaviside step function. The factor | ^ |r^ ; 1 1 accounts for the flux of approach- 
ing needles (therefore requiring the factor 6 (— ^|r^|))- The following two step functions 
and the delta functions enforce the conditions |f > \sij\, |f > |s_jj|, and |r^| = + for a 
collision, assuming no lateral extension of the needles. Finally fry, acting on a phase space 
function A(T), replaces linear and angular momenta before a collision by those after collision 
according to the rules given in Eq. One can equally consider the evolution backwards 

in time. The corresponding Liouville operator is given in Ref. |32 . 



IV. PROJECTION OPERATOR FORMALISM 



In the following, autocorrelation functions of phase space functions or observables A(T, t) 
will be calculated: 

$(0 = (A*(T,0)A(T,t)) 

= (A*(T J 0)exp(iC + t)A(T,0)} (28) 

The bracket (...) defines an average with respect to the canonical ensemble with weight 
oc exp(— Ti/T). A* denotes the complex conjugate of A. For later purposes, the Laplace 
transform is defined as follows: 

r°° 1 

*(z)=i/ dte*4(t) = -{A%r,0)——A(r,0)). (29) 
Jo z + L + 

The function $(z) is analytic for Im(z) > 0. 

Following Mori and Zwanzig we can represent the autocorrelation function in terms of 
a continued fraction, which is suited to an approximation by truncation. Let us denote 
the resolvent in Eq. ()29|) by 1Z = (£+ + z)" 1 and define a scalar product of two dynamic 
variables A and B via the equilibrium expectation value as (A\B) := (A*B). The Laplace 
transform of the autocorrelation can then be written as a matrix element of the resolvent 
*(z) = -(A\(z + C+)- 1 A). 

The idea behind the Mori Zwanzig formalism is to consider explicitly the dynamic evo- 
lution in a subspace of dynamic variables and treat the dynamics of the remaining variables 
in a simple approximation. The crux of the matter is to choose the right subspace. Popular 
candidates are the conserved quantities, which are special because of their long relaxation 
times in the hydrodynamic limit. The fast degrees of freedom are then approximated by 
a single relaxation time, similar in spirit to the single collision time approximation of the 
Boltzmann equation. We shall in the following consider the subspace spanned by p(r,u), 
but first review the general formalism. 

We consider a set of dynamic variables {A{\ with normalisation matrix (xo)ij — 
(Ai\Aj) = (A*Aj). The projector onto the subspace spanned by the {Ai} is given by 
P = ^2 i j\Ai)(xo 1 )ij(Aj\. The orthogonal projector is denoted by Q = 1 — P . Fol- 
lowing the Mori-Zwanzig formalism we obtain an equation for the correlation functions 
<&ij(z) = —(Ai\(z + C + )~ l Aj). In matrix notation it reads 

(z + QoXo 1 + M (z)xo 1 ) Hz) = -xo (30) 
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with (Oo)y = (Ai\C + Aj) and the memory kernel 

(Mohiz) = -(A i \£ + Q (z + Q £ + Q )- 1 Q £+A j ). (31) 

Since the memory kernel is again a matrix element of a modified resolvent taken with the 
state Bi = Q C + Ai, we can iterate the procedure, thereby representing the autocorrelation 
as a continued fraction. We define a second projector f\ = \B i )(xi 1 )ij(Bj | and Q\ = 
1 — Pi and apply the same procedure to M (z), 

(z + fiar 1 + M 1 (z)xl 1 ) M (z) = -xi, (32) 

with (xi)ij — (Bi\Bj), (^i) = (Bi\C + Bj) and a second memory kernel M^z). 

In principle this procedure can be iterated an arbitrary number of times. However the 
computation of the static expectation values Xi an d becomes increasingly complicated, so 
that one ususally truncates the continued fraction after one or two steps. In the case of hard 
core interactions the "restoring forces" f2j are in general imaginary, because the Liouville 
operator is non hermitean. Hence even a truncation of the continued fraction with M\ = 
- the approximation scheme that will be adopted in this paper - will give rise to damping 
effects. This is in contrast to differentiable potentials with a Hermitean Liouville operator. 

V. CORRELATIONS OF COUPLED DENSITY ORIENTATION FLUCTUA- 
TIONS 

A. Observable and matrix elements 

We are interested here in the anisotropic motion of a tagged rod-like macromolecule. In 
particular we want to compute the anisotropic diffusion constants D± and D\\ as well as the 
rotational diffusion constant Dr (see Eq.(£Q)) from a kinetic theory approach. The dynamic 
variables of interest are thus the density of a tagged particle p(r, u) = 5(r — r )5(u — u ) for 
all orientations u. It is convenient to consider its Fourier transform p(k, u) = e~ lkr °S(u — Uq) 
and compute the matrix of correlation functions 

Su.»'(M = -(p(k,u)|(z + £ + )- 1 p(k, U / )) (33) 

where the orientation vector u is the "matrix subscript" . Actually u is continuous, so that 
the matrix equations turn into integral equations. To apply the formalism of Mori and 
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TABLE I: Geometrical factors cm, c± and cr for homogeneous rods (W = 6) and "dumbbells" 
(W = 2), rounded to 3 decimal places. 





smooth rods 


rougl 


l rods 




W = 6 


W = 2 


W = 6 


W = 2 


c ll 








10.631 


10.202 


C_L 


7.253 


8.626 


11.802 


13.390 


CR 


0.544 


0.684 


0.871 


1.058 



Zwanzig, we identify 

A u = p(k, u) = e- 4kro 5(u - u ). (34) 

We consider two steps in the continued fraction and only discuss the simplest approximations 
Mi = 0. As will be shown below, the resulting correlation function $ uu /(k, z) has the right 
hydrodynamic limit. It furthermore accounts for free particle motion approximately. More 
elaborate approximations as well as the effects of futher steps in the continued fraction will 
be discussed in the conclusions. 

The computation of susceptibilities Xi an d "restoring forces" ( 1 is straightforward but 
cumbersome and hence delegated to appendix A (see also Refs. The results of this 

calculation are 



(Xo)u,u' = {A U \A U ,) = (4tt) ^(u - u 
(floW = (A u \£+A u> ) = 
(Xi)u', u = (A u \A ul ) 

Hill \ 

n 



L(L k i + Li 2 ) s(u - u') 

4tt U I 



(35) 
(36) 

(37) 



(fii) U ',u = (A u \£ + A u ,) 



rp \ 3/2 



5/2 \ m J 



C ||L 2 (k • u) 2 + c ± L 2 {k 2 - (k ■ u) 2 ) + 4W 2 c R L 2 5{u' - u) 

(38) 



Here the angular momentum operator is defined as usual (see Eq. ©). The constants cm, 
cj_, and cr are geometrical factors specific to needles, depending on the dimensionless ratio 
W = mL 2 /(2I). They are defined in App. [B]and their numerical values are listed in Tab. U 
for W = 6, corresponding to homogeneous rods, and W = 2, corresponding to rods whose 
mass is concentrated at the endpoints ("dumbbells"). 
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B. The correlation function 



Given the matrix elements obtained above, we proceed to find a solution for $ u ,u' (k, z) 
within the approximation M\ = 0. Since Qq and Xo are trivial quantities, we will concentrate 
on Qi and Xi an d drop the subscript 1 from now on for notational simplicity. Additionally, 
we absorb a factor of 1/xo = 4-7T by defining ip = 4tt§. Multiplying Eq. (|30|) from the left 
by (z + fix -1 ) (arid using Eq. ([32))). we are then left with the matrix equation (w.r.t. u and 
u') 

{z 2 + zQx" 1 ~ 4vr X V = -(* + ^X~ l ) (39) 

What is an appropriate set of functions to represent <p in? The matrix x is diagonal in 
the eigenfunctions of the isotropic Laplacian, i.e. the spherical harmonics. The matrix Q is 
diagonal in the eigenfunctions of the anisotropic Laplacian, denoted by yi m {ic, u) and related 
to the oblate spheroidal polynomials (or oblate spheroidal wave functions) Si m (ic,r)) in the 
same way as the spherical harmonics are related to the associated Legendre functions, 

yi m (ic, u) = yi m (ic, r), 4>) = -===S lm (ic, r])e tmcl> , (40) 

v 2tt 



with rj = u ■ h/k. The properties of spheroidal wave functions are discussed in App. O 

As will become evident below, we do not need to keep the dependence on the azimuthal 
angles <fi and in Eq. ()3*9*j) and integrate the equation over these variables to obtain 

J drf' {z 2 + zQx' 1 ~ 47rx) w »</VV = ~(* + ^X~ l )r,rf- ( 41 ) 
Here <p m > is defined as 

Vrm' = 7^ J ' d<p J d(j)'<p uu > (42) 

and analogously 

fi W = £liSi (ic, T))Si (ic, rj') 
i 

Xvv' = ^XiSio(0,v)Sio(0,v')- 
i 

The expansion coefficients are given by 



(43) 



n ' = m (w^ + A ' o( " c2) 

1 f T T 
47r V m 1 



(44) 
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with c* = ^(L*)' and. = ^(£) 3/2 . 

In App.[0]we show how to diagonalize y> m > perturbatively for k — > 0. For the time being, 
however, we choose the set S/o(0, 77) oc -P^ 7 /) to expand 



^^^o(0,r/)^ (0,r/'). (45) 



/./' 



In particular, we will later need the matrix elements <poo an d <f2o, which are calculated from 
the diagonalized tp in App. [DJ The result is 

^00 = 5- ztt^ : I 46 ) 

„ c 2 -n 2 

V?20 = 47T 



9^ (z 2 + 2X0 ^0 - 4ttxo) (z 2 + % - 47%) 
^02 + 0(k 4 ), (47) 



where only terms up to and including 0(k 2 ) have been retained. 



C. Numerical simulations 



We have performed simulations of the system of hard rods using an event-driven algorithm 
as described in Ref . I35I (with no inelasticity). We used system sizes of N = 800 to N = 2048 
and densities between p = 0.25 and po = 100- F° r the most part we used homogeneous 
rods with W = 6. The results of the simulations will be discussed below along with the 
theoretical ones. 



VI. RESULTS 



In this section we are going to discuss the results of the approximate kinetic theory 
obtained by setting Mi = 0. One might expect that our theory reduces to the diffusion 
equation of Perrin in the limit of small wave numbers kL <C 1 and long times. However, 
such a straightforward limit does not exist as we will show below. 

There are two special cases where the comparison between microscopic theory and the 
Perrin equation is straightforward. One case concerns purely rotational motion. We show 
that the theory predicts a crossover as a function of density from weakly damped free rota- 
tions to rotational diffusion. We furthermore obtain explicit expressions for the rotational 
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diffusion constant of the Perrin theory as well as for the damping of the free rotations. The 
other important case which allows for direct comparison with the Perrin theory is purely 
translational motion, where we recover the isotropic diffusion constant. 

The third and most interesting case, however, concerns the coupling of translational and 
rotational degrees of freedom. The comparison with the Perrin diffusion equation is not 
directly evident. In fact, we derive - within kinetic theory - the full time dependence of 
the mean square displacements of a tagged rod parallel and perpendicular to its initial 
orientation. We then compare our expressions to the corresponding ones obtained from the 
Perrin theory. An important problem arising in the analysis is the appearance of short 
time scales present in the kinetic theory which are not easily decoupled from the long- 
time behavior. One would like to map the latter to the Perrin diffusion equation which 
involves only long-time diffusive time scales. Even within the Perrin theory, the anisotropic 
translational motion characterized by the diffusion constants D» and D±_ is only perceptible 
on small times (compared to the rotational diffusion time) before the anisotropy is hidden 
by the long-time isotropic translational diffusion of the center of mass of the tagged rod. 
On the side of kinetic theory, the time interval for anisotropic diffusion to be observed is 
bounded from below by ballistic motion which is isotropic a priori. We comment on this 
point in more detail below when comparing our analytical results with simulations. 



A. Rotational diffusion 

We first discuss purely rotational motion, as obtained from the results of Sec. |V] by 
setting k = 0. In order to compare with simulation results later on we focus on the Fourier 
transform tp(uj) instead of the Laplace transform <p(z) of the correlation function. Since in 
the time domain tp(t) is a real symmetric function, the Fourier and Laplace transforms are 
related via <f(co) = 2lm(<p(z = uj + z0)) where to is real. 

The correlation function for k = is diagonal in the angular momentum quantum num- 
bers and is given by 

m ' \z 2 -ujf\ 2 + \z\ 2 C 2 + 2lm((z*(z 2 -uf)) 1 ' 

with uf — + 1) and £ = Ann-!-. The Fourier transform is 



= 2 T-2 ifc ( 49 ) 
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The correlation function Eq. (J49|) has two poles in the lower complex frequency plane 
at u = — «C/2 ± a/ (uf — C 2 /4) and the corresponding complex conjugates in the upper 
half plane. It shows qualitatively different behaviour depending on the reduced density 
Po = uqL 3 . For low density the poles have a nonzero real part, corresponding to weakly 
damped oscillatory motion of the rods: 

fuit) = — exp {-(t/2) (u t cos (uit) + C sin (u^t)) (50) 



The real part of the poles vanishes above a critical density p cr it(0 — [}(]> + 1) Jiy 



1/2 



which depends on /. For large density the time dependent correlation function is a sum 
of two exponentials with the long relaxation time Ti ong (/) = (DrI(1 + 1)) _1 and the short 
relaxation time r s hort = C^ 1 - The decay on long times corresponds to rotational diffusion 
and is in agreement with the Perrin theory. The rotational diffusion constant is given by 

Dr = ^J L (51) 
PoLc R V m 

Our approximate theory does not describe the free particle limit po ~~ * correctly. In 
particular the weakly damped oscillations acquire additional damping due to free particle 
motion (sometimes called motional narrowing). We discuss a simple way to account for 
these additional terms approximately. 

If there are no interactions, the angular correlations are simply given by 

r°° j ji z 
ifl(t) = (Pi(cosut)) th = duju- Pi(cosujt) e~~. (52) 

Jo 1 

For example for / = 1 (the case we will need below) the correlation function ^\i{z) can be 
expressed in terms of the exponential integral function Ei(z)^. 

A simple way to incorporate free particle motion is the following. One represents the 
correlation function ipu as a continued fraction with Q = but M/0 (see Eq. (|30|)). 

z+Mi/xi 

Then the memory kernel is expressed in terms of (p^, 

m°(z) 

Z+M ^' = 4 * X <^MT1' (54) 



1(3 



and can now substituted in the continued fraction for the interacting system. This procedure 
yields 

§M + i) + <M>x, 
^ = -^M + u + W (55 » 

We now compare the Fourier transform (pu(u>) = 2lm(ipu(z = uj + iO)) for the case 
/ = 1, which corresponds to the correlation function <pn(t) = (u (t) • u (0)), with numerical 
simulations in Fig El Agreement is very good for small densities because the free particle 
limit has been incorporated. Our approximate theory also reproduces the crossover from 
damped oscillations, corresponding to a peak at finite frequency in the spectrum, to purely 
relaxational motion, characterized by a peak at zero frequency. 

B. Translational dynamics 

Next we discuss purely translational motion, obtained from the general results by setting 
/ = /' = 0, i.e. by examining the correlation function (p 00 as given in Eq. ()46|) . Its Fourier 
transform is given by 

^ ooH = \^-{kv^y + e^ (56) 

the thermal velocity v% h = T/m. 



with ^ = v ■ Here we have introduced the Enskog collision frequency v — i^y^f- and 



The corresponding time dependent correlation function is the sum of two exponentials. 
The decay rate of the fast component, tq, is determined by the collision rate according to 
Tq 1 = u(cn + 2c_i_)/(27r) 2 . The decay rate of the slow component is diffusive as it should 
be, due to particle number conservation. Since we consider translational motion only and 
consequently have set I = in the equation of motion, we should recover the isotropic 
diffusion constant of the Perrin theory. The latter predicts D lso = (D\\+2D±)/3. Comparison 
with the Perrin theory allows us to identify 



, th (2vr) 2 Ttt 6nL 



As ° ~ v C|| + 2c ± ~ V m p (c|| + 2c ± ) ' (5?) 

The ratio of the rotational to the isotropic translational diffusion constant is given in the 
Kirkwood-Riseman theory^ by Dr/D iso = 9/L 2 . Even though this theory is designed for 
a rodlike molecule in a solvent where hydrodynamic interactions are important, we do find 
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almost the same result for a melt of smooth needles with a homogeneous mass distribution: 
Dr/D iso = (c|| + 2c_i_) / (3crL 2 ) = 8.89/L 2 (see Tab. |TJ). The ratio does depend on the mass 
distribution though and is different for smooth and rough rods. 

In summary, our approximate kinetic theory correctly describes particle diffusion which 
dominates the long time decay of the translational correlation function. It yields an ap- 
proximate expression for the isotropic diffusion constant, which is inversely proportional to 
the density po- The decay of the nonconserved variables has been approximated by a single 
relaxation time, which is inversely proportional to the Enskog collision frequency. 

C. Coupled translational and rotational dynamics: kinetic theory 

1. From correlation functions to mean square displacements 

Our next goal is to calculate the mean-square displacement of the tagged rod, decomposed 
into the directions parallel and perpendicular to its initial orientation u (0). We define the 
tensor 

(Ar a (t)Ar^(i)Mo 7 (0)^(0)), (58) 

where Ar(t) = r (t) — r (0) is the displacement of the tagged rod. This tensor R^s is 
isotropic and symmetric under exchange of the indices a and (3 and also of 7 and 5. These 
symmetries imply that there are only two independent tensor components Ri and R2 such 
that 

Ra/3-yS = SaP^&Rl + (SaeySpS + S a sSfS 1 )R 2 . (59) 

In order to calculate Ra/3-yS, we express it in terms of the correlation function $, 



rj2 

Raw = -QkQk- <e Jk - Ar( V(0)u 05 (0)) 



k=0 







(60) 



dk a dk 



ft 



dudu u w 5 $ uu /(k, t) 



k=0 



The total mean-square displacement can easily be expressed in terms of this tensor as 

3 

, RaaW = UXl + 6/^2 = - 2^ 7^2 / du d\l'<P uu , (k, t ) ^ = -3 ° 

a/3 



o2 r 

(Ar 2 (t)> = R ^(3/3 = m + 6R 2 = - J du du '®™' ( k ' *) 



- o; 2 ^00 

k=o ok z 



k=0 

(61) 
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(recall that $ and </? are identical up to a factor of Aw). In the last step the derivatives 
become simpler since y?oo is only a function of k 2 (cf. Eq. f!46|) ) . 

In terms of R\ and _R 2 the parallel mean-square displacement is given by (cf. Eq. (|59|)) 

(ArJ(t)) = £ i? a/3Q/3 = 3Ri + 12i2 2 . (62) 

In order to evaluate this further, we need to know one more equation besides Eq. (jfiljl to 
determine R\ and R 2 . This is provided by 

^™ = 3jR i + 6i? 2 = - du rfu '«) 2$ ™' ( k ' *) I k= 

d 2 , 2 I 

'^5^02 + </%)) 



9k 2 y/5 I k=0 

Note that each term in the sum over a does not depend on a because of symmetry and the 
fact that ip 00 and y?02 only depend on k 2 . Thus it suffices to evaluate the term with e.g. 
a = 3. 

Eqs. fEQ), (E2J, and dS3) then result in 

d 2 

(M(0) = faoo + V^ra) , (64) 



<% 2 

a 2 



fc=0 



(Arl(t)) = -^(2^00 - ^02) fe=Q . (65) 

Anisotropic diffusion from kinetic theory 

Given the expression derived above for the mean square displacement parallel and perpen- 
dicular to the initial orientation, the last step remaining is the calculation of the correlation 
functions (p Q0 (k,t) and <P2o{k,t) in the time domain from the corresponding expressions de- 
rived in the previous section for the Laplace transform (depending on z). Therefore, we need 
to analyse the poles of v 3 oo(^, z ) and <f2o(k, z). The denominators of ^00(^5 z) and <p2o(k, z) 
consist of factors of the form z 2 + —z — iiiXn (n = 0, 2). The two roots of each of these 
expressions are at 

Zn,± = -7^(1 ± VTTlG^M). (66) 

While <poo(k, z) only contains the n = factor and thus has the two poles z ,±> l fi2o(k, z ) has 
both the n = and n = 2 factors and therefore has all four poles. 
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Given these poles, we can expand y?oo an d y?2o from Eqs. (j4T)j) and (f4T)) into partial 
fractions and carry out the inverse Laplace transforms. From Eqs. (164}) and (J63j) we then 
get the mean-square displacements. Before we show the results of this straighforward but 
tedious calculation, we note that the long-time behavior on long length scales (k — > 0) 
is determined by the pole or poles closest to the origin. Each pole z n ± corresponds to a 
timescale l/(iz nt ±) (provided z n> ± is purely imaginary, see below). Inserting the values from 
Eq. ()44j) these timescales are (immediately setting k = where possible) 

1 m fm 6Ln 3 / 2 

iz , + T V T Po(2c ± + C||) 

1 [7n27i 3 / 2 L 



T > = ^rH^-J 1+ ^- iSnW ° w4)] (68) 

1 /^2tt 3 / 2 L 



^ " = ir^wT R I, 1 " V 1 -**>M*<%>) (69) 

i z o,~ D iso k 2 

Generically, only one of these times is long, namely r 3 , which is proportional to 1/k 2 . For 
large densities, however, T2 becomes proportional to the density and thus large, too. This 
indicates that, at least for high densities, there are two long time scales. The other two 
times scales, To and t±, are short. The first one, r , is related to the Enskog collision time. 
The second one, T\, is of the same order as r . For large densities it agrees with the short 
rotational timescale r s h or t from Sec. IVI Al It characterizes the crossover from free rotational 
motion to rotational diffusion. 

After performing the inverse Laplace transforms, we find 

(Ar 2 ) = 6Aso (t + r (exp(-t/r ) - 1)) (71) 
(Ar 2 ) = i(Ar 2 ) + (Ar 2 niso ) (72) 



3 

(Ari) = ^(Ar 2 )-(Ar 2 niso ) (73) 



where 



47Tfc(c ± -C||)L 2 v 

{ anW " zw 2 c R -^n^v^-i/v) 



(74) 



fc=0 

Due to the density dependence of Q 2 the poles z 2 ,± acquire a real part for small densities, 
such that the pole is no longer purely diffusive but oscillatory. Note, however, that the mean- 
square displacements are real quantities whether z 2) ± are purely imaginary (corresponding 
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to an exponential decay of the nonlinear part of the mean-square displacements) or not 
(corresponding to an oscillatory behavior). 

The anisotropic mean square displacements calculated from kinetic theory have the fol- 
lowing important characteristics: 

• For short times t C ro one obtains ballistic behavior 2(Arjj(£)) ~ (Arj_(i)) ~ 2t> 2 h t 2 
The time r is identical to the fast timescale calculated in Sec. IV1 Bl 

• For long times £ ^> r 2 the mean square diplacements approach the isotropic center of 
mass diffusion behavior 2 (Ar (?(£)) ~ (Ar^(t)) ~ AD iso t 



For large densities, the pole Z2 - is given by 

" - 1 - 8t I) = - udr (75) 

with the rotational diffusion constant Dr given in Eq. (JoTj) . The timescale r 2 is 
therefore r 2 = (6-D^) -1 = Ti on g(^ = 2) where ri on g is the long rotational timescale 
(see Sec. IV! A|) . It also agrees with the crossover time from Perrin theory, see Eqs. 
and fllDj) . 

For intermediate times, when Tq <^ t ^ t 2 (provided such an interval exists, which is 
only the case for high enough density), there is a linear (in t) regime in (Ar 2 niso ) which 
corresponds to anisotropic diffusion. The diffusion constants D± = D[ SO — \AD and 
Dn = D iso + AD in this regime can be calculated from Eq. (|74j) . resulting in 

2ttk(c± — cii)L 2 



AD 



3W 2 c /? (l/r 2 -l/r 1 )(l/r 2 -l/r ) 
D R 



(76) 

po^oo (C_L - C\\)L 



(2c ± + c {l )W 
3. Anisotropic diffusion from simulations 

The results of our analytic calculation are compared with the simulations in Figs. El ID 
and 03 for homogeneous rods (W = 6) and in Fig. [U] for dumbbell molecules (W = 2). 

The agreement between theory and simulations is very good for small densities po < 10 
but at these densities there is hardly any discernible anisotropy (see Fig. OJ). 

For higher densities an anisotropy window opens up between the parallel and perpendic- 
ular mean-square displacement curves. This is exemplified in Fig. 0] For densities above 
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w 20 the isotropic diffusion constants (slopes of the mean square displacements) increase 
with density in the simulation. This effect is well know n 23 1 24 and is due to enhanced parallel 
diffusion along a tube formed by surrounding particles. This involves correlated particle 
collisions and is thus not captured within our theory. Therefore the theory does not fit the 
data quantitatively but Fig. 0] shows that the anisotropy is very well captured qualitatively. 
But even though anisotropy is clearly present, it can (at these densities) not adequately be 
described by anisotropic diffusion constants since the anisotropy window is still small and 
within it, no convincing "straight line fit" appears possible. This is true for both theory and 
simulation. 

Fig. shows a direct comparison of theory and simulations of the mean-square displace- 
ments for densities po — 40 and po — 100. Here it becomes obvious that the isotropic 
diffusion constant is underestimated by the theory. The simulations at po = 100 now very 
clearly show anisotropic diffusion with a well-defined perpendicular diffusion constant D±. 
The parallel diffusion constant Dh, on the other hand, is still not well defined even at this 
density. 

The anisotropy becomes even more pronounced for dumbbell molecules with W = 2 at 
p = 100, which is shown in Fig. The reason the effect is more pronounced for dumbbells 
lies in the H^-dependence of the diffusion constants. While AD is proportional to 1/W (see 
Eq. JZOJO, Mso is only weakly dependent on W through c\\ and cj_. Therefore anisotropic 
diffusion is more pronounced for smaller W. 

VII. DISCUSSION AND OUTLOOK 

Before we go into a discussion of the details of our results we will give a brief overview of 
our main findings. The aim of this paper is to gain some microscopic understanding of the 
transport properties of a gas of thin hard rods beyond the phenomenological Perrin equation. 
Purely rotational and translational dynamics of a tagged particle can be described by simple 
rotational or translational diffusion, in agreement with the Perrin equation, allowing us to 
calculate the rotational and isotropic translational diffusion constants. The motion of a 
tagged particle along and perpendicular to its original orientation shows marked anisotropics. 
The nature of the anisotropy is, however, very different from the one predicted by Perrin 
theory, at least for low and medium densities. This is confirmed by simulations. The 
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Perrin equation predicts anisotropies which can be described as anisotropic diffusion for all 
densities. In contrast, our results for low densities exhibit only a very drawn-out crossover 
from ballistic motion for short times to isotropic diffusion at long times with no intermediate 
regime which could be characterized as anisotropic diffusion. For medium densities po ~ 40 
such a regime begins to emerge; it is however still far from being well-defined. Even for 
densities po = 100 where the mean-square displacement is clearly anisotropic, anisotropic 
diffusion is only seen in the perpendicular motion. 

We have derived an approximate theory for the coupled translational and rotational 
motion of thin hard rods, starting from a microscopic model. The dynamic evolution is de- 
termined by binary collisions, in which the normal component of the velocity of the contact 
point is reversed and the tangential component is either left unchanged (smooth rods) or 
reversed as well (rough rods). The Mori-Zwanzig projection operator formalism has been 
applied to the correlation of a tagged rod's density p(r, u, t), specifying its position r and 
orientation u at a given time t. Matrix elements of the collision operator have been com- 
puted, resulting in a generalized Enskog theory, in which only uncorrelated binary collisions 
are taken into account - similar in spirit to a "StoBzahl Ansatz" of the corresponding Boltz- 
mann equation. Two steps in a continued fraction expansion give rise to an integral equation 
for the correlation function. This equation has been solved in the limit of small wavenumber, 
where it is rendered diagonal by linear combinations of the spherical harmonics. We have 
also performed numerical simulations, using an event driven algorithm for hard needles. 

Purely rotational motion (k = 0) is as expected and in agreement with previous discus- 
sions. For small density we observe damped oscillations with a frequency u 2 = T/I. In 
the dilute limit it is essential to include free particle motion, which provides the dominant 
damping mechanism. 

For large densities the dynamics is purely relaxational. The crossover happens at a critical 
density p„it(0 = 87r 3 /(/ + 1)/(crW) which depends on the angular momentum I, the moment 
of inertia via W = ^^f- and on whether we are considering smooth or rough rods. For 
Po > Pcrit(0 our approximate kinetic theory predicts two relaxation times. The longer one is 
determined by the rotational diffusion constant according to T\ ong (l) = (DrI(1 + 1)) _1 with 
the rotational diffusion constant _D# given by Eq. (J51|) . It depends on the mass distribution 
along the rod and is different for smooth and rough needles. The rotational diffusion constant 
Dr is always larger for smooth rods than for rough rods, e.g. for a homogeneous mass 
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distribution along the rod we find f) s ™ oth j £)™ ugh « 1.60. 

The long time relaxation T\ ong (l) = (/^/(Z+l))" 1 is in agreement with the Perrin equation, 
where Dr enters as a phenomenological parameter. In addition our theory also predicts a 
fast decay with a relaxation rate r s hort = Dpi IT or alternatively IstmX = -%^/(7 + l) = (^ k ) 2 . 

''"long J- ^PO 

This short relaxation time is mainly determined by the moment of inertia and hence roughly 
three times larger for dumbbells than for a homogeneous mass distribution. 

Next we discussed purely translational motion (1 = 0). We find relaxation with a short 
and a long relaxation time, the latter being diffusive 77 = (Disok 2 )^ 1 and in agreement with 
the Perrin equation, which allowed us to identify the isotropic diffusion constant D lso as 
in Eq. (|57|) . It depends on the mass distribution and is different for smooth and rough 
needles. The isotropic diffusion constant is always larger for smooth needles than for rough 
ones, e.g. for a homogeneous mass distribution along the rod we find Df™ oth /D r °^ « 2.36. 
Comparing different mass distributions we always find a larger diffusion constant for the 
homogeneous mass distribution than for the dumbbells. In addition to the diffusive long 
time decay, the kinetic theory predicts partial decay on a microscopic timescale, which is 
given by r = D iso m/T. 

The most interesting case is of course the coupled translational and rotational motion, as 
described by the time dependent correlation $ u , u '(k, t). This function has been diagonalised 
in the long wavelength limit by linear combinations of spherical harmonics, which for general 
k, l,m differ from the solution to the Perrin equation. We attribute this difference to the 
fact that the Perrin equation is not simply a gradient expansion but involves the angular 
momentum operator, which does not contain any small parameter allowing for a systematic 
expansion. 

Anisotropic diffusion is best discussed with the help of the particle's displacement vec- 
tor, projected onto or perpendicular to the particle's initial orientation, i.e. (Arjj(t)) and 
(Ar\(t)). In general, we find three dynamic regimes. For short times t <C To the motion is 
ballistic. For the longest timescales, i.e. times larger than the timescale for rotation of a 
needle, t ^> 1/(6-Dr), isotropic diffusion prevails. In between a time regime appears which 
is characterized by anisotropic diffusion: a linear increase of the mean square displacements 
(Ar|(i)) ~ D\\t and (Ar\(t)) ~ D±t with anisotropic diffusion constants D\\ and D±. To 
clearly see this intermediate regime, the density has to be sufficiently large, such that the 
timescale for rotational diffusion of the rods is long. In the simulations these three time 
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regimes are clearly visble for densities po > 40 (see Fig. HJ). The analytical theory correctly 
predicts anisotropic diffusion qualitatively (see Fig. Hj), however for such high densities no 
quantitative agreement between simulations and theory is achieved, see Fig. |3J Anisotropic 
diffusion is more pronounced for a higher moment of inertia, the highest value being reached 
for dumbbells with all the mass concentrated at the endpoints of the needle (see Fig. |BJ). 
Lowering the density leads to shrinking of the intermediate time regime, such that for small 
densities (po < 10) one only observes a crossover from ballistic motion to isotropic diffusion. 
For these densities the agreement between theory and simulation is very good, as shown in 
Fig. El 

Since the intermediate anisotropic regime is only present for high densities, we conclude 
that the Perrin equation is an adequate description of the dynamics of thin hard rods only 
for high density. In order to derive it from a microscopic theory, it is not sufficient to let 
the wavenumber become small; in addition the rotational diffusion constant has to be small 
or in other words the rotational timescale has to be long. 

A perspective for future work based on the microscopic theory presented here would be to 
exploit the correlation function to calculate intermediate scattering functions. These could 
then be measured in depolarized light scattering experiments. Moreover, the microscopic 
theory is valid for all wave vectors. However, up to now, a diagonalization procedure used 
to derive explicit expressions for tpui exists only for small wave vectors. 

APPENDIX A: MATRIX ELEMENTS 

In this appendix, the essential steps necessary to calculate the matrix elements xo, 
Q (u',u), Xi(u',u), and fi 1 (u',u) are presented. 

1. Generating functional 

The structure of the collision operator requires a specific coordinate system for the trans- 
lational momenta (or velocities). 

In fact, the operator iC, appearing in matrix elements is reduced to the iV-fold operator 
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for collision of needle with needle 1, 

i£ + = N\V- uj 0(-V • u ± a)6 - \s 01 \\ 6 - \s 10 \\ 5 (|r^| - + ) (&+ - l) . (Al) 
The equality 

— |r^- 2 | = u ± • V r .sgn(r 12 ■ u±) (A2) 

has been used, and a is short for sgn(ri2 • uj_). The relative velocity of contact points V 

and its representation in the orthonormal basis Uj_, Uq, uq have been defined in section ITTT1 

Before performing the transformation, let us recall the definition of the integration measure 

in average (...) in terms of canonical coordinates: 

i 

dT = Y\_dr i d(p i d9 i dp ri dp ipi dpg i (A3) 

i=0 

The normalization reads as 

Z-i = J_J: 1 L_ (AA) 

V 2 (4vr) 2 (2vrmT) 3 (2tt/T) 2 1 ; 

Next, the canonical coordinates are transformed to new ones which are standard normally 
distributed (with respect to the Boltzmann factor exp(— flH.)). Changing moreover transla- 
tion momenta p rj to CMS momentum and relative momentum, the transformation reads as 
follows: 

* = 72k (p °- pi) (A5) 

7 = 7^ (P0 + Pl) (A6) 

Pe t = -/Lp, 8 (A7) 



P Vi = 1_ V % . (A8) 



with % — 0, 1. Positions are transformed according to 

r i = r - ri (A9) 
r = r (A10) 

The resulting integration with respect to roi will be split into the following three components: 

c = r i • u ± (All) 

» = '->= 7== § ( A12 > 

y/l - (U ■ Ux) 2 

. = -s m = r 01 ■ u - (u °; U '» (r "' U ° ) (A13) 

y/1 - (U • Ux) 2 
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This decomposition complements the one for the relative velocity of contact points. The 
corresponding Jacobian is given by 



The relative velocity of contact points given in terms of the new coordinates reads as follows: 



As the velocity V r and relative position r i have been decomposed with respect to the 
base (uj_, Uq", uo), the expression given just above will be rewritten likewise. Following a 
suggestion by Huthmann, one may transform from standard normally distributed variables 
(P^PflJ with respect to the base (e^e^) to new standard normally distributed variables 
(vi,Wi) with respect to the base (u^,Uj_), giving the following expression for V r : 




(A14) 




(A15) 




(A16) 



Later, on will be eliminated using 



u 



l 



— Au + Bu { 



i 



(A17) 



where 




(A18) 




(A19) 



So a decomposition of V r with respect to the base (uj_, u^~, u ) results. 

In terms of the variables (7, Xi u o, wo, v%, w\) the 2-needle Hamiltonian reads 



pn = - ( 7 2 + x 2 + «o + ^0 + vl + w{) . 



(A20) 



Accordingly, the integration measure multiplied by the normalization changes to 



1 



1 



1 



drdadbdcd{l dQid~fdxdv dw dvidw 



(A21) 



Z~ l dT 




1 
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where dQi = dipisin(6i)d8i defines the integration with respect to spatial angles. Now, the 
collision operator involves a projection of V r . Moreover, 7 is related to V r by Eq. (|A16j) . 
Therefore, 7 will have to be expressed in terms of V r in the integration measure and the 
Hamiltonian. The integration measure involving a factor by dx = (^) 3 ^ 2 dxdydz, the 
Hamiltonian is 



m 
2T 



m 



(x 2 + y 2 + z 2 ) + (v 2 + w 2 ) (l + a 2 ^) 



?) (l + b 2 



w 



m\ m 
2l) +X T 



-y—^ljiBbvx + avo) 



m 



(awo + bwi 



j(Ab Vl ) 



m ~ 



m 



+ ab—Bv vi + ab—woWi + 7 



(A22) 



The variable x will require a particular treatment as will be shown later on. Moreover, the 
variable 7, denoting the CMS momentum, will not appear in matrix elements involving only 
relative translational momenta, and therefore will be treated seperately. Defining the vector 
Xby 

X T = (2/,js,Vo,u>o,Vi,u>i) 



the Hamiltonian may be rewritten in the following form 



-X T AX 




(aw + bwi] 



;1 



(A23) 



(A24) 



where A is the following matrix: 
/ 

A = 





m 
2T 





2T\ I 





^sJ^fbB 


\ 







m 
2T 













m 
2T 


Ft 
\ l a 





(l + a 2 §) 





abf T B 
















(l + a 2 f) 





ab% 


m 

2T\ 






abf T B 





(l + 6 2 §) 
















ab§ 





(l + 6 2 §) ) 



(A25) 



Let us now consider the generating functional with respect to 

3 j h Vl , h Wl ) , 



(A26) 
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namely 

G(H) = J dXexp (-PH + H T X) 



(A27) 



The integration is gaussian and the result is: 



IT / m \ -V 2 

G(H) = (27T) 3 — (1 + -(a 2 + 6 2 ) ) exp (-^ (A28) 
m \ 21 J 



with the effective Hamiltonian 



^ cff = Y (If) (1 + ^2 + 62)) " ^A^H + 1 7 2 
2T(l + §(a 2 + o 2 )) V / 1 lJ 



The inverse matrix A 1 is given by 

(l + i(« 2 + ^ 2 )) -f(i£6 2 ) -y^a -y^oE 

^(ABb 2 ) ZL(i + %pA 2 ) jzM 

fa 10 



(A29) 



o oo = , m ;. 2 2 f^ o 



2/ 



f&5 x/t^ 6 10 

afef l+a 2 f 
x U l+f(a 2 +6 2 ) U l+|i(a 2 +6 2 ) / 

(A30) 

Whenever matrix elements involving the collision operator appear, the delta function 
S (l r oil -0 + ) =5 (|c| - 0+) in terms of new position coordinates as well as the dependance 
on sgn(c) need to be taken care of. In this case the following relation is useful: 

/oo rco 
dc5 (|c| - 0+) .../(... , sgn(c)) = lim / dc (S(c + e) + 6(c - e)) .../(... , sgn(c)) 
-oo e— »0 J „ OG 

•••(/(•••• -1) • ./•(••••-!)) 

(A31) 
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2. Matrix elements 



Matrix elements explicitly involving the collision operator can be easily calculated when 
represented as moments of the variables (x, y, z, v , w , v±, W\, 7) using the generating func- 
tional derived above. In fact moments involving (y, z, v , w , vi, Wi) can be obtained by 
replacing expression involving a variable say y to some power by the respective derivative 
with the order equaling the power. The integrations with respect to x and 7 are treated 
separately. 

In fact, a matrix element (A*iC' + A) usually has the following form 

f L/2 f L/2 

(A*iCA) =C da db 

J-L/2 J-L/2 

x I ... I JiQoMe-W 
x (6(-x) + O(x)) |x| 

x /(fi , Qi,x, y, z, v , w , v 1 ,wi, 7) (A32) 

where 

C = ^VN(m/2T) 3/2 (A33) 

Zj 

The function /(Ho, tt\,x, y, z, Vq, wq, v±, w±, 7) contains the functional form of the factors A* 
and A in the matrix element after (b^ — l) has been carried out. Other factors arising from 
the collision operator have been written separately. The factor C is a normalization factor 
which results after integrating with respect to the CMS position r and transforming form x 
to (x, y, z). The generating functional now permits reducing the number of integrations as 
follows: 



(A*iCA) = C [ ' da f ' db I ... [ J(fi ,fii)(e(-x) + 6(x))|x| 

J-L/2 J-L/2 J Jn ,Ui,x,~f 

x f(n , Q u x, 7, d hy , , d hz , d hvo , d hwo , d hvi , d hwi )G?(H) (A34) 

In our case, no dependances on 7 exist, so the corresponding integration may be immediately 
carried out, giving a factor of (27r) 3 / 2 . 

Let us outline the calculation for the matrix element ^(u',!!). It involves 2 matrix 
elements (B*CBi) for i = 1,2 resulting from the two terms in A u . Taking into account the 
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change of Alio an d Ap after a collision (which are given for the general case by Eq.s (J18|) 
and (EH)), the matrix element reads as follows: 



'•-L/2 f* f* 

{B*CB 1 ) = (-i)C / da db ... J 



X 



e m [v. • i ^V r + a\^j{v Q v^ + w u ± ) + b^j^u^ + u>iu_l) j 



x |x| (0(-ar) + 0(x)) — k • (aui + 7iu + 72^) S(u' - u )<5(u - u ) (A35) 



and 



-L/2 i-L/2 

(B* 2 CB 2 ) = (-z)V^V^C I da db I ... I J 



L/2 J —L/2 J J Qo,fl\,x,y,z,VQ,WQ,vi,w 



am _ Fjl _ ( I , T IT < T . < 

x — e pri I u (u • V r ) + u (u ■ u x )6W yt>i - V r - aW j(^ou + w u_l) - OW y (vii^ + wi 

x |x| (9(-x) + 6>(x)) (aui + 7 2 u L ) ;y <5(u' - u )<5(u - u ) (A36) 

The coefficients a, 71, and 72 are defined in Eq.s (fTTj) and (fl§j). and contain dependancies 
on x, y, z, (uo • ui), as well as a and b. 

Using the generating functional as above one finds the following results, 

1 /j-\3/2 

(BJXfli) = I- L 2 [ C || (k • u) 2 + c ± (A; 2 - (k • u) 2 )] 6(u' - u) (A37) 

(B* 2 CB 2 ) = m -i- (-) ^ ^ 4 f ^) c R L 2 5(W - u) (A38) 



8tt 5 / 2 VW W 2 

The coefficients cm, c±, and cr contain the remaining integrations over a, b, and the relative 
angle between Qq and fii, and are defined in appendix B. Adding the two matrix elements 
gives the expression for fii(u', u) in Eq. (|3HI) . 

APPENDIX B: SOME GEOMETRICAL INTEGRALS 



Let us first define 

f l/2 

da I db (Bl) 



1/2 f l/2 fll 



1/2 



1/2 v/l + W(a 2 + b 2 ) 
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The constants appearing in Eq. (j38|) are given explicitly as follows where W = mL 2 /(2I) 

■1/2 pl/2 



2tt / da f db (1 + W(a 2 + b 2 )) 1 ' 2 

V 2 / 7-1/2 J -1/2 

x r dx (1 - x 2 Y' 2 (l + Wa 2 + Wb 2 x 2 ) 

J.! 1 J (l + VT(a 2 + & 2 )+WW(l-X 2 )) ^ 



= 7T 

•1 



Jo + f ^1 7T / da / do (l + W(a 2 + o 2 )) : 

V ^ / .7-1/2 .7-1/2 



,2,1/2 (1 + ^(1-X 2 )) 

l l J (l + W{a 2 + b 2 ) + W 2 a 2 b 2 {\-X 2 )) 



Cfl = vr 2 / 2 + P-^) 7T / da/ doa 2 (l + W/(a 2 + o 2 )) 

V 2 / y_i/ 2 7-1/2 
x C d ri - 2,1/2 (i + ^ 2 (i-x 2 )) 

7_! 1 XJ (l + jy(a 2 + 6 2 ) + W/ 2 a 2 o 2 (l-2 2 )) 
They have been evaluated numerically for W = 6 and W = 2 (see Table H} 



(B4) 



APPENDIX C: SOME PROPERTIES OF SPHEROIDAL WAVE FUNCTIONS 



For a general introduction to spheroidal wave functions we refer the reader to Refs. 1291130 . 
In this work we only need spheroidal functions of either prolate or oblate type. The prolate 
functions are defined as the solutions of the eigenvalue equation 



d " 2, d , X m2 



+ WO - " CV) s,UC ,) = o, (Ci) 

where / > and m are integers obeying —l<m<l and £ is a real constant. The oblate 
functions can be obtained from Si m ((,r]) by replacing ( by i(. We prefer working with an 
orthonormal basis and therefore choose the normalization 

J dr,St n (C,v)Sv m (C,v) = Su' (C2) 



which is unfortunately different from the one used in Refs. |2£ 

The Si m ((,r]) can be expressed in terms of associated Legendre functions P\ m as 



Nim(0 

Here dJ, m (C) are expansion coefficients and 



SlUC,V) = -^Yl^OPr+rnAv)- (C3) 



^o^/E'(w 2r+2 2 m+1 (r+ r r )! ( C4) 
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is the normalization constant. For methods to obtain \i m (() arid d l ™ from Eq. ()C1|) we again 
refer the reader to Refs- I^lls^ . For our purposes it is sufficient to know the eigenvalues A/ m (C) 
up to order ( 2 and the scalar product 

(£WC),<W0)) = J drjS; m (C,v)Si' m (0:V) (C5) 

up to order ( 4 . The results of this tedious calculation (which is omitted here for the sake of 
brevity) are 



a^(C) = a / { 2 + c 2 aS + ^(C 4 ) (C6) 

!(/ + !) -2 
(21 + 3) (21 



'('+i)+c 2 2 '<l+l-!'" 2 1 :W ) (C7) 



and 



(Si m (C), S Vm {0)) = 6 W + Cg m iv + C%mi> + 0(( 6 ) (Gl 



where 



yf(l - 1 - m)(l - 1 + m)(l - m)(l + m) 

9mW = — 



(C9) 



(CIO) 



2(21 - l)V(2/-3)(2/ + l) 
^/(Z' - 1 - m)(/' - 1 + m)(Z' - m)(Z' + m) 
l+2 ' 1 ' 2(2/'- l) 2 v^2/'-3)(2/' + l] 

7 1 - 5w , . s 

A; m — A;, m 
A (2) - A ( , 2) 1 

-4^ TWdmlV - 2$ U' (9ml, 1+2 + 9ml,l-2) and 

A; m — A Fm 

_ ~ 1 + m)(l - 1 - m)(f + m)(l - m) 
(21- 1V(2Z+ l)(2/-3) 

Here it is implied that all coefficients are if their indices lie outside their range of validity 
1,1' > and —I, —I' < m < I, I'. For some values of the indices the coefficients appear to 
have vanishing denominators. However, in these cases the prefactors are 0, too, and the 
respective terms should be taken to be 0. 



APPENDIX D: DIAGONALIZATION OF $ 

In order to diagonalize the matrix of correlation functions if vri >(h, z), it is sufficient to di- 
agonalize B m > = z(x~ 1 ) m ' + (x ^X -1 )^'- Since Xrm' * s diagonal in the basis of (normalized) 
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Legendre polynomials y I + \Pi{tj) = Si (0, rj) and VL m i is diagonal in the basis of Si (ic, rj), 
where c 2 = (c± — c\\)L 2 k 2 /(4W 2 cr) is small in the limit k 2 — > 0, we aim at a perturbative 
solution. The starting point is the representation of B in the basis {S^O, 77)}. The matrix 
elements are 

B lv = J drjdrj'S lo (0,rj)B w S llo (0,rj'). (Dl) 

Since x is already diagonal in this basis, nothing needs to be done for the term z\~ x and 
we get z{x~ X )w = zSw/xi- 

The term x" 1 ^" 1 is more complicated. Inserting 1 = J2i \Sio(ic))(Sio(ic)\ in the two 
places between x _1 and Q we get 

B lv = z8 lv / X i + {x^X~ X )iv (D2) 

= zS w I xi + S w ^ + AB} V + AB 2 , + AB? V , where (D3) 
Xi 

ABj, = c 2 g m ,^ l - D4 

XiXi' 

AB 2 V = —(h wo —5 lo + h m —5 VQ ), (D5) 
Xo Xr Xi 

c 4 

A-B«/ = gojigoji'ttj (D6) 

X/X*' ^ 

to order k 2 . The calculation is straightforward but care has to be taken with xo an d ^0 
as they are of order k 2 while all other eigenvalues of x an d ^ are °f order 1. Therefore 
the matrix elements AZ?q 2 and AB\ Q are of order 1. However, all other off-diagonal matrix 
elements of B are of order k 2 or smaller. The diagonal element B 00 is of order 1/k 2 while 
all other diagonal elements are of order 1. In order to be able to use perturbation theory we 
employ a Jacobi transformation 

B B 

Rw = 5 a > + (cos — l)(<W^'o + Si 2 Si> 2 ) + sin — {5 m 5 V2 - fe<W- (D7) 

which makes RBR T = D + k 2 A where D is a diagonal matrix and A is a perturbation. Note 
that R is really only required up to order k 2 which is 

Rw = Sa> + ^(SioS V2 - 5 l2 5 m ) + 0(k A ). (D8) 

-DOO 

The eigenvalues bi of B (up to order k 2 ) are given by the diagonal entries of RBR T since 
the Jacobi transformation merely changes the basis vectors but not the eigenvalues. The 
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diagonal entries of RBR T are equal to b\ = Bu = z / xi + ^i/xf for all / 7^ 0, 2. For the two 
exceptions one gets 60 = -Boo + 2B 2 2 / B 00 (/ = 0) and 62 = -B22 — B 2 2 /B 00 (/ = 2). These 
are, however, higher order corrections which are irrelevant for our purposes and we may 
therefore use b\ = z/xi + ^i/xf f° r ai l I- 

Finally, the eigenf unctions S'i (77) can be obtained through perturbation theory. The 
matrix Bu> is symmetric but not normal (i.e. [B, B'] 7^ 0, except if z happens to be purely 
imaginary). Therefore B and of course also RBR T do not have an orthogonal system of 
eigenvectors and we have to calculate both left and right eigenvectors perturbatively. This 
procedure yields for the matrix O r of right eigenvectors of RBR T (0] v is the Z'th component 
of the /th eigenvector), 



0,W„. + (1-M^>£, (D9) 

Ol — Of 

while the matrix O 1 of left eigenvectors is 

o\ v = 5 IV + (1 - su>) f™^ = my (Dio) 

and the asterisk denotes the complex conjugate. 

The right and left eigenf unctions £[(77) and Sj(r)) of the original operator B are then 

Sf(v) =^(0^)^^0(0,77) (Dll) 
v 

S 1 l (r ] ) = J2(O 1 R)u>S l/0 (0, V )- (D12) 
v 

The product 1 R is, up to order k 2 , 

(O l R) w = 8 IV + ^(5 l0 8 V2 - S l2 S V0 )+ (D13) 

-DOO 

As an application which is needed in the main text we calculate some matrix elements of 
the correlation function ip in the basis of spherical functions 5/0(0, 77). We have 



- ^(^ (0)|^)-^(^|^o(0)) (D15) 

m bm 

- J2(O r *R) ml —^(O l R) ml ,, (D16) 



z h 

On 
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and from this and Eq. (|D14|) we get 



1 



z + Qq/xo 



(D17) 




z(z + fio/xo) - 47TXo' 



omitting terms of higher order than k 2 . Likewise we obtain 




B02 1 , B 02 1 



(D18) 



-n 2 



(D19) 



9v/5 (z 2 + zxo^o ~ 4ttxo) {z 2 + % - 4vrx 2 ) ' 
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FIG. 1: The collision plane Eqi spanned by the orientations of two rods. 

FIG. 2: Orientational correlation function <pn(uj) from theory (Eq. (|55|l) and simulation for various 
densities po- 

FIG. 3: Mean square displacements and for homogeneous rods (W = 6) at densities po = 4 
and po = 10 from theory and simulations. The uppermost line indicates the exact ballistic behavior 
for small times. Distances are measured in dimensionless units with L = 1. 

FIG. 4: Mean square displacements and for homogeneous rods (W = 6) at densities po = 40 
(simulation) and po = 100 (theory). At these densities, an anisotropy window opens up between 
the parallel and perpendicular mean-square displacements. 

FIG. 5: Mean square displacements and r\ for homogeneous rods (W = 6) at densities po = 40 
and po = 100. 

FIG. 6: Mean square displacements for dumbbell molecules with W = 2. 
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